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ABSTRACT 

We suggest novel statistics for the CMB maps that are sensitive to 
non-Gaussian features. These statistics are natural generalizations of the 
geometrical and topological methods that have been already used in cosmology 
such as the cumulative distribution function and genus. We compute the 
distribution functions of the Partial Minkowski Functionals for the excursion set 
above or bellow a constant temperature threshold. Minkowski Functionals are 
additive and are translationally and rotationally invariant. Thus, they can be 
used for patchy and/or incomplete coverage. The technique is highly efficient 
computationally (it requires only O(N) operations, where N is the number of 
pixels per one threshold level). Further, the procedure makes it possible to split 
large data sets into smaller subsets. The full advantage of these statistics can be 
obtained only on very large data sets. We apply it to the 4-year DMR COBE 
data corrected for the Galaxy contamination as an illustration of the technique. 

Subject headings: cosmic microwave background, cosmology, theory, 
observations. 
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1. Introduction 

Observations of the Cosmic Microwave Background (CMB) provide valuable 
information about the early Universe. In addition to the possibility of measuring the 
cosmological parameters the CMB data can provide very important constraints on the type 
of the seeds that led to the structure formation (see e.g. Pond fc Jaffe 1998p . Inflationary 
theories predict Gaussian density perturbations with nearly scale-invariant spectrum (e.g. 
turner 1997] and references therein). The Gaussianity of the density perturbations results 



in the Gaussianity of the CMB temperature fluctuations at the surface of last scattering. 
Thus, testing the Gaussianity of the CMB fluctuations becomes a crucial probe of inflation. 
On the other hand, if the seeds for the structure formation are due to topological defects, 
such as strings and textures (e.g. Prandenberger 1998| and references therein) then the 



non- Gaussianity of the temperature fluctuations may probe fundamental physics at high 
energies. 

However, even if the temperature fluctuations were Gaussian at the surface of the last 
scattering they may acquire small non- Gaussianity due to subsequent weak gravitational 
lensing ( see e.g. |Seljak 19961 , |Bernardeau 19971 , IWinitzki 1998|) as well as due to various 



astrophysical foregrounds (see e.g. Pandy et al 1996 ). Higher resolution maps (MAP 



PLANCK) will make it even more problematic. 

Establishing the Gaussian nature of the signal is also important for practical reasons: 
Some current techniques for estimating the power spectrum are optimized for the Gaussian 
fields only (e.g. [beldman et al 199~^| , see also the discussion in |Knox et al 1998| and |l'erreira 



et al 1998| ). 



Standard tests for non-Gaussianity are the three-point correlation function or 
bispectrum and higher order moments. However, in practice negative results of the 
non-Gaussian tests can hardly be conclusive since only infinite number of n-point correlation 
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functions can prove that a field is Gaussian. A distribution may appear to be Gaussian up 
to very high moment and then be non-Gaussian ( [Kendall fc Stuart 1977|) . The high-order 



correlation functions are very expensive computationally for the large data sets (0(N m ), 
where m is the order of the correlation function and N is the number of pixels). Thus, any 
statistic that is sensitive to non-Gaussianity and computationally efficient is very useful. 

Here are other well known examples of the tests sensitive to non-Gaussianity: peak 
statistics (|Bond fc Efstathiou 1987] , |Vittorio fc Juskiewicz 19871 , [Novikov fc Jorgensen| 



1996|) , genus curve (integral geometric characteristics) QMelott et al 1989| , |Coles 1988 



INaselsky & JNovikov 1995Q ; global Minkowski Functionals (hereafter MFs) ( |Gott et al TM) 



Schmalzing fc Gorski 199"B , Winitzki fc Kosowsky 19UH ). These functionals also have been 



considered for CMB polarization field by |Naselsky fc Novikov 1998] . Minkowski Functionals 



( Minkowski 1903]) were 'properly' - i.e. in the context of differential and integral geometry 



introduced into cosmology by |Mecke, Buchert, fc Wagner 1994] as a three-dimensional 



statistics for pointwise distributions in the universe and then for the isodensity contours of 
a continuous random field by [Schmarzing fc Buchert 1997] . We will discuss the MFs in the 



following section. Here we just mention that in the two-dimensional case of the temperature 
maps the global Minkowski Functionals are the total area of excursion regions enclosed 
by the isotemperature contours, total contour length, and the genusQ or the number of 
isolated high-temperature regions minus the number of isolated low-temperature regions. 
Partial Minkowski Functionals are the same quantities but used as characteristics of a single 
excursion region. 

Kogut et al. (1996) measured the 2-point and 3-point correlation functions and the 
genus of temperature maxima and minima in the COBE DMR 4-year sky maps. They 
concluded that all statistics were in excellent agreement with the hypothesis of Gaussianity. 



Tn flat space the genus is equal to the Euler-Poincare characteristic. 
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Uolley et al 1996| measured the genus of the temperature fluctuations in the COBE DMR 



4- year sky maps and came to a similar conclusion. Heavens (1998) computed the bispectrum 
of the 4- year COBE datasets and concluded that there was no evidence for non-Gaussian 
behavior. 

However, Ferreira, Magueijo & Gorski (1998) studied the distribution of an estimator 
for the normalized bispectrum and concluded that the Gaussianity is ruled out at the 
confidence level at least of 99%. 

The first analysis of two-dimensional theoretical maps of the temperature fluctuations 
that used the total area, length of the boundary and genus for the excursion set was done 
by |Gott et al 1990| although without referring to Minkowski Functionals. Then |Schmalzmg 



fc Gorski 1998 discussed the application of the MFs to the COBE maps stressing the 
importance of taking into account the curvature of the celestial sphere (the manifold 
supporting the random field of the temperature fluctuations). [Schmalzing fc Gorski 1998 



actually applied the statistics to the COBE data and argued for its advantage as a test of 
non-Gaussian signal. They concluded that the field is consistent with a Gaussian random 
field on degree scale. 



Salmi, Sathyaprakash & Shandarin 1998 suggested that the Partial Minkowski 
Functionals (PMF) may be used as quantitative descriptors of the geometrical properties of 
the elements of the large-scale structure (superclusters and voids of galaxies). In particular, 
they argued that two simple functions of the PMFs called the shapefinders can distinguish 
and reasonably quantify the structures like filaments, ribbons and pancakes. 

Here, we describe a new statistical tool which is sensitive to non-Gaussian behavior 
of random fields. We suggest using the distribution functions of the Partial Minkowski 
Functionals and the number of maxima for a given threshold level. As an illustration we 
apply it to the full sky temperature maps obtained by subtracting the galaxy contributions 
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from 4- year COBE observations ( [Bennet et al 1992 , Bennet et al. 1994 ) The main goal of 



the paper is the demonstration of the potentiality of this novel technique. 

The outline of the paper is as follows. In Section 2 we briefly review the Minkowski 
Functionals. In Sec. 3 we outline the numerical algorithm used for the evaluation of the 
Minkowski Functionals. Then, as an example of application this method we present the 
first results of analysis of the 4-year COBE maps. Finally, in Sec. 4 we discuss the results 
and potentiality of the method. 



2. Minkowski Functionals 

Let us consider one connected region Ri of the excursion set with 
v{6,ip) = (AT(9,{p)/T)/(To > u t , where o"o = ((AT/T) 2 ) 1//2 , u(9,(p) is the measure 
and v t is the threshold. If the region is complex then it may require very many parameters 
to fully characterize it. However, we consider only three parameters: the area of the region, 
a,, the length of its contour, Zj, and the number of holes in it n^. These are three partial 
Minkowski Functionals. In the following analysis we also count the number of maxima of 
u(9,(p). To obtain the global Minkowski Functionals, we compute these quantities for all 
disjoint regions of the excursion set, i.e. taking the sums A = Ea^, L = EZj and G = Hgf 
"number of isolated v > v t regions" — "number of isolated v < v t regions". The last 
quantity is the genus introduced into cosmology long ago by |Doroshkevich 1970 and pott et 



|al 1990| . The total area A(u t ) is clearly proportional to the cumulative distribution function 
of the random field. At high (low) threshold levels the excursion regions appear as isolated 
hot (cold) spots. 

The Minkowski Functionals have several mathematical properties that make them 
special among other geometrical quantities. They are translationally and rotationally 
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invariant, additive Q, and have simple and intuitive geometrical meanings. In addition, it 
was shown flHadwiger 1957|) that all global morphological properties (satisfying motional 



invariance and additivity) of any pattern in D- dimensional space can by fully characterized 
by D + 1 Minkowski Functionals. 

Global MFs of Gaussian fields are known analytically; in the two-dimensional flat space 
they are: 

L(u) = ±ue-4, (1) 

..2 

(2tt) 3 / 2 20j ''" 



G{y) = TOTok^e v 2 . 



X 

where = J e~ x 'dx' is the error function. Their dependence on the spectrum can 
be expressed only in terms of the length scale of the field 8 C = ^ where ctq and G\ can be 
calculated from the spectrum C\ 

<7„ = l/47r£i(2J + l)C„ 

<T 1 = l/4nZi(2l + l){l + l)lC l . 

The analytic formulae for the partial Minkowski Functionals are not known even for 
Gaussian fields. However, it does not impose a principal obstacle in their application since 
they can be calculated numerically. It is worth stressing that in practical application, in 
addition to the mean value of some quantity one has to know its variance. In most cases 
the variance is not available in an analytic form even if the mean value is. For instance, 
one can calculate analytically the mean number of hot/cold spots but the variance of this 
number can be estimated only numerically. 



2 In particular, additivity means that the Minkowski Functionals of the union of several 
disjoint regions can be easily obtained if the Minkowski Functionals of every region is known. 
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3. Application to Two-dimensional Maps 

We identify all disjoint regions above a given threshold v > v t for positive peaks and 
below the threshold v < — Vt for negative peaks. For every region i?, we compute three 
Minkowski Functionals: 

The area, v\ = 

The perimeter v\ = U i.e. the length of the boundary 
The number of holes - equivalent to genus - v\= gi 
The number of maxima v\ = n m i within the region. 

We then study the cumulative distribution functions F(u t ,v k ) (k = 1,2,3,4) of these 
quantities. 

3.1. Data 

We used the DMR 4-year whole sky maps, where all Galactic emission was removed. 
Two independent methods were employed to separate the Galactic foreground from the 
cosmic signal. The maps' construction is described in detail in two papers (Pennet et 
al 1992j , Pennet et al. 1994j) , they were released in the DMR Analyzed Science Data 
Sets (ASDS). The two techniques are: One method is the so-called combination method 
(map 1). Here they cancel the Galactic emission by making a linear combination of all 
DMR maps, then cancel the free-free emission assuming a free-free spectral index and 
finally normalize the cosmic signal in TD temperature. The subtraction method (map 2) 
constructs synchrotron and dust emission maps and subtracts them from the DMR data 
sets. The Galactic free-free is then removed. The analysis presented in this paper was done 
for both maps and the results of the analysis were similar. 
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3.2. Numerical Algorithm 

In this subsection, we describe the numerical algorithm for calculation of the 
distribution of the partial Minkowski functionals on the sphere and application of this 
algorithm to the COBE data, considered in the previous section. 



3.2.1. Maps simulations. 

In our simulations we use spherical coordinate system to assign pixels on a sphere. 
Here we consider the temperature distribution on the pixelized map as the function of two 
variables in the coordinate system: — it/ 2 < 9 < n/2 and — ir < ip < it. In fact, this function 
is defined only in the points (9^, y?fc 2 ), so that: 

Vk u k 2 = v(0ki,<Pha), 0*i = hhe, Vk 2 = hK ■ ( 3 ) 

We also assume, that he = h v = h = || where M is the number of pixels in the ip direction. 
The total number of pixels is, therefore, M 2 /2. 

The original COBE maps have been recalculated according this pixelization in the 
following way: 

AT data (6, <p) = Bj AT cobe (6', <p')e ^dcos(9')d<p' (4) 

where ATqobe and AT data are temperatures defined in the points of COBE cube pixels and 
in the points defined by the variables in Eq. (3) above, 7 is the angle between the pixels, 
70 = 7° is the smoothing angle and B is the normalization. 

The temperature fluctuations are completely characterized by the spectrum coefficients 
C™. Using this description one can write the following expression for the temperature of 
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the relic radiation: 

oo m=l 

AT data (0,v) = E E Carrie,?), 

l=2m=-l (5) 

VdataxyiV) — Tat 5 — 

* data' 

where YJ m are the spherical harmonics. The summation in Eq.(5) is from 1=2. The 
term with 1=1 is the dipole component. This term has been removed from the COBE 
data before analysis because the contribution of this term in the AT fluctuations can 
not be separated from the contribution due to the motion of the observer relative to the 
background radiation. 

We have simulated 1000 different Gaussian realizations of the temperature distributions 
on a sphere so that we can compare the distribution of partial Minkowski functionals in the 
observational data with a random Gaussian field. This was accomplished in the following 

way: 

oo m=l l 

!=2m=-( (Q) 

where a™ are independent random Gaussian numbers with zero mean ie. (a™) = and with 
unit variances ((a[™) 2 ) = 1. The power spectrum C\ in (6) were obtained using Cf 1 from (5) 
as follows: 

m=l 

Y. (Cf'f 



3.2.2. Calculation of Partial Minkowski Functionals 

Let us introduce some threshold v t on a pixelized map of the CMB. Each isolated 
hot-spot (regions with v > u t ) can be considered as the cluster characterized by the 
area, boundary length and Euler characteristic or equivalently by genus (both are directly 
related to the number of disjoint boundaries). For example, the total area of the map 
where v > v t is the sum of the areas of all isolated hot-spots, where v > v t . The global 
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Minkowski Functionals, i.e. the total area, total boundary length and total genus can be 
found by summation of their partial values over all clusters on the map. Computing partial 
Minkowski functionals on the pixelized map we require that the algorithm satisfies the 
following convergence properties: 

v l k \ p — v\— > 0, as h^O 

and (8) 
(fcl)^O(^) A; = 1,2 

where v\\ p denotes k-th Minkowski functional of i-th cluster, calculated on the pixelized 
map and v \ is the exact value of this functional on the continuous field. In the following 
analysis we use linear interpolation, so that z — 1 for our algorithm. 

Pixels (&4, k 2 ) inside the regions where v > u t satisfy the condition v^^ 2 > u t . We define 
pixel (ki, k 2 ) inside this region as the inner boundary pixel if the field is below the threshold 
level v t at least in one of its four neighbors {{k\ + 1, k 2 ), (k\ — 1, k 2 ), k 2 + 1), (k± : k 2 — 1)) 
(eg. fkx+iM < v t) see Fig- 4. We approximate a smooth boundary curve by the polygon 
using linear interpolation of the field between inner and outer boundary pixels (Fig. 4) and 
thus finding the intersection of the boundary curve with the grid lines: 

e b = hh + h — UklM , Vb = k 2 h (9) 

^1+1^2 — u ki,k 2 

for if grid lines and 

6 b = hh, ip b = k 2 h + h — UklM — (10) 

^1^2+1 — v kiM 

for 9 grid lines. 9 b and ip b denote coordinates of the boundary points X m = (9 b , ipb) on the 
polygon. This polygon obviously converge to the smooth boundary line as h — > 0. 

Then, the algorithm of cluster analysis consists of three steps. 

• 1) Identifying the boundaries and computing their lengths. 

First, we search for closed boundary lines of the level v = u t . Then, each set of 
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— * 

boundary points is ordered by letting X^ +1 to be the nearest boundary point to the 

— * 

point X£. The length of a closed boundary line is: 

m=M n +l 
m=l 

where M n is the total number of boundary points in the n-th closed line of threshold 
level (X^ n+1 = X™) and the norm 

\X m+1 - X m \ = [(9 m+1 - 6 m f + S in 2 (^±i±^)(0 m+1 - 0J 2 ] 1 / 2 . 

The first point X\ is arbitrary. Different boundary lines in the map correspond to the 
arrays of boundary points (X™) and inner boundary pixels (Y£). The total boundary 
of an isolated region v > v t may consist of a number of closed lines (two lines in 
Fig.4). 

• 2) Finding all the boundaries of a connected region (cluster), computing the total 
boundary length and genus. 

We combine all closed lines which are the boundaries of the same cluster by using 
arrays of inner boundary pixels (Y£). Suppose, we wish to check whether two different 
lines are the boundaries of the same cluster or not. These lines correspond to two sets 
of inner boundary pixels Y^ 1 and Y^ 2 . If we take two arbitrary inner pixels one from 
each set and connect them by a path along grid lines (see Fig. 4) then the path can 
intersect the boundaries N\ nt times (i = 1,2), where N\ nt > 0. If the both numbers 
Nl nt and Nf nt are even then both inner boundary pixels belong to the same cluster 
otherwise they belong to two different regions (clusters). Therefore all boundary lines 
which belong to one cluster form its boundary with the total boundary length equal 
to the sum of their lengths. The number of closed lines for each cluster is equivalent 
to the genus of this cluster. Thus we find total number of clusters and two partial 
Minkowski functionals for each of them - the length and genus. 
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• 3) Computation of the areas of clusters. 

All pixels situated between the inner boundary pixels of a cluster belong to the 
cluster. The area of the cluster can be roughly approximated by the total area of all 
these pixels (including inner boundary pixels). 

3.3. Results 

Figures 5 and 6 show the cumulative distribution functions F(u t , Vk) (where 14 = J2i v l) 
for the two COBE maps described above along with the mean Gaussian values and variances. 
The mean and variance were obtained from 1000 random realizations of Gaussian fields 
having the same amplitudes as shown but different sets of random phases. In both Fig. 5 
and 6 we see significant deviations from Gaussianity. It may be interesting to note that 
each statistic shows the greatest difference with the Gaussian value at different thresholds: 
F(A) at v t = -0.5, F(L) at u t = -l, F(G) at v t = ±1, and F(N m ) at v t = 0, 0.5, -1. These 
differences are roughly the same for both maps and suggest that each of the four statistics 
carry different statistical information. 

As one might expect, more detailed information can be obtained from the partial 
Minkowski Functionals. Figures 7 through 11 show the partial Minkowski Functionals at 
ten thresholds v t = ±2, ±1.5, ±1, ±0.5 and ±00. Every figure shows two curves, one for 
the positive {v > v t , solid lines) and one for the negative [y < u t , dashed lines) thresholds, 
having same absolute magnitude v t for each map. Thick and thin lines correspond to COBE 
maps 1 and 2 respectively. The mean Gaussian curve (that obviously does not depend on 
the sign of the threshold) is the dotted line and la Gaussian variance is shown as a shaded 

3 The thresholds v t = ±0 and —0 correspond to the excursion sets v > and v < 
respectively. 
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area. 

The main features of Fig. 7-11 are the following: 

• Figure 7, v t = 2: The functions F(a) and F(l) show strong non-Gaussian signal; F(g) 
and F{n m ) roughly consistent with the Gaussianity. 

• Figure 8, v t = 1.5: All statistics shows strong non-Gaussianity. 

• Figure 9, v t = 1.: The strongest non-Gaussian signal comes from the distribution of 
maxima, F(n m ); other statistics are roughly compatible with Gaussianity. 

• Figure 10, v t = 0.5: All statistics show marginal disagreement with Gaussianity. 

• Figure 11, v t = 0: All statistics are in rough agreement with Gaussianity. 



4. Discussion 



We suggest new statistics to test the Gaussianity of CMB maps. These statistics 
are the distribution functions of the Partial Minkowski Functionals of the excursion set 
for a given threshold level v t . We also compute the distribution function of the number 
of maxima in the isolated region. The Partial Minkowski Functionals have transparent 
geometrical and topological meanings of : 1) the area, 2) the perimeter and 3)the genus of 
each disjoint region. Introducing these statistics is a natural step further. It generalizes the 
statistical techniques used before: the cumulative distribution function that correspond to 
the total area of the excursion sets, global genus ( poroshkevich 1970| , |Gott et al 1990Q , the 
total length of the boundary (|Gott et al 1990|) . The set of three characteristics mentioned 
above is known as global Minkowski Functional and has already been used in cosmology 
flMecke, Buchert, fc Wagner 19941 , ftchmalzing fe Buchert 19971 , |Winitzki fc Kosowsky"T9"9"8" 



|Schmalzing fc Gorski 1998|) . |Sahni, Sathyaprakash fc Shandarin 1998| argued that the 
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Minkowski Functional can describe the morphological characteristics of isolated regions or 
objects. 

It is well known that Minkowski Functionals are invariant under translations and 
rotations and are also additive. Due to these properties the MFs can be used in cases 
of incomplete or patchy coverage. The additivity allows to analyze large data set by 
splitting them into a set of smaller subsets. Measuring the partial MF is very efficient 
computationally: for a fixed threshold it requires only O(N) operations, where N is the 
number of pixels. 

Clearly, computing the distribution functions of the partial MF require much larger 
data sets than we used here for an illustrative purpose. However, the upcoming MAP and 
Plank missions will provide the necessary resolution. Needless to say that by knowing the 
partial MFs one can easily obtain global MFs. Thus, these statistics completely incorporate 
the global MFs. 

We show that the reconstructions of the whole sky maps obtained by subtraction the 
Galaxy contributions from the COBE maps are strongly non-Gaussian. This was suspected 
by many authors who used only the polar caps for the analysis of the cosmological signal 
(e.g. |Colley et al 1996] , |Schmalzing fc Gorski 1998| , |Ferreira et al 19981) . It is possible that 



some of the non- Gaussian signal in our analysis is due to errors in the Galaxy removal, 
since the Galactic signal is obviously non-Gaussian. For instance, |Colley et al 19961 found 



no deviations from Gaussianity in the genus curve while we see significant non-Gaussianity 
in our genus measurements. Unfortunately, by analyzing the Galaxy caps only we use 
only about half the data and in addition split the remaining half into two disjoint regions 
affected by the boundaries. Since the largest number of disjoint regions is about ten at 
v t = 1.5, 1 (Fig. 9,10) in the whole sky, the analysis of the caps only does not seem to be 
reasonable. This paper provides a feasibility study in the strengths and usefulness of the 
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Partial Minkowski Functionals, not a definitive results regarding the non-Gaussian nature 
of the COBE-DMR maps. 
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Fig. 1. — The COBE map 1 constructed using the combination method (see §3.1). The 
heavy lines show 0, la and 2a (the interior of 2a contours is shown in black); the light lines 
show —la and —2a (the interior of — 2a contours is shaded). 
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Fig. 2. — The COBE map 2 constructed using the subtraction method (see §3.1). The 
notations are same as in Fig. 1 



-22 - 




Fig. 3. — An example of a Gaussian map with the same as in COBE map 1 amplitudes 

1(1 + 1)C,/(T 
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Fig. 4. — A region bounded by the level contours v > u t is shown. The closed polygons are 
the approximations to the boundaries based on the linear interpolation. Stars and circles 
are two sets of the inner boundary pixels corresponding to the boundaries. The dashed 
line shows a possible path on the grid connecting a pair of the inner boundary pixels that 
may belong to the same region or two different regions. If such a path intersects the both 
boundaries (i.e. the boundaries corresponding to the pixels in question) an even number 
of times then both pixels belong to the same (connected) region, otherwise they belong to 
different regions. In this case both numbers are zeros (i.e. even) and thus they belong to 
the same region. 
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Fig. 5. — The cumulative distribution functions of global Minkowski Functionals and the 
numbers of maxima/minima as a function of the temperature threshold (given in units of a 
for COBE map 1. Circles and triangles show the values for positive and negative thresholds 
respectively. The error bars corresponds one a dispersions calculated in 1000 Gaussian 
realizations. 
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Fig. 6. — Same as Fig. 5 but for the COBE map 2. 
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Fig. 7. — The cumulative distribution functions of partial Minkowski Functionals for both 
COBE maps. The shaded strips show ±a regions for Gaussian realizations. Solid line shows 
the cumulative distribution functions for positive (y > is t ) and dashed line for negative 
{y < i> t ) thresholds. Thick lines correspond to COBE map 1 and thin lines to COBE map 
2. The threshold v t = 2a. 
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Fig. 8. — Same as Fig. 7 with v t = 1.5a. 
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Fig. 9. — Same as Fig .7 with v t = la. 
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Fig. 10. — Same as Fig. 7 with v t = 0.5a. 
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Fig. 11. — Same as Fig. 7 with v t — 0. 



